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I. INTRODUCTION 

A major source of the thermal pollution of rivers, lakes 
and estuaries is the hot water discharges from power plants. 

As the demand for electric power in this country continues to 
increase, more and more power plants will be built. The 
trend is that large power plants, including more nuclear 
powered stations, will be built, and groups of these units 
will be located at a single site - usually in coastal areas. 

This will result in sharp increases in thermal pollution. 

Thus, proper management of the water resources and the related 
thermal pollution problems is essential to preserve the en- 
vironment in estuarine and coastal marine waters. In order 
to achieve this goal it is necessary to have a thorough under- 
standing of the motion and diffusion of thermal discharges 
and the extent of the region affected. To obtain the under- 
standing of the hydrodynamics and thermodynamics of earth's 
water basins requires extensive field data collection. Even 
the most rudimentary data acquisition field trip is extremely 
expensive. Experimental observations from physical models 
are often unsatisfactory owing to the impossibility of obtain- 
ing appropriate similarity. Developing theoretical predictive 
ability is considered to be the better way to obtain the large 
volume of information required to understand the environment. 

Models can be used predictive ly to arrive at meaningful 
decision parameters a priori for site selection for power 
plants and other modifying factors. 

In a NASA sponsored feasibility study (refs.l and 2) completed in 
January 1974, a team of researchers at the University of Miami 
concluded that it is feasible to develop a generalized three- 
dimensional predictive, mathematical model, involving remote 
sensing and in-situ measurements to supply necessary para- 
meters. It was recommended that the development of the mathe- 
matical model be conducted over four phases covering a time 
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period of three years. Three sites were proposed to be used 
for the development and testing of the model. These three 
power plant sites are the Florida Power and Light Company's 
Turkey Point and Cutler Ridge facili ties, and Port St. Lucie 
Nuclear Power Plant units on Hutchinson Island. These three 
sites are located in South Florida, close to both the Kennedy 
Space Center and the University of Miami. Another interesting 
feature is their geographical contrast. The Turkey Point and 
Cutler Ridge facilities are located at a shallow lagoon type 
estuary, while the Hutchinson plant discharges into the off- 
shore continental shelf. These sites should provide a good 
test area for the generalized nature of the proposed three- 
dimensional model. The initial plan was that a single model 
would be developed for application to thermal pollution prob- 
lems caused by thermal discharges with a six-month Phase I 
study to conduct preliminary model development. A twelve- 
month Phase II study was to be devoted to model developing 
using Biscayne Bay site as a test case. A twelve-month Phase 
in study was to apply the model to the Hutchinson Island site. 
A six-month final phase study was to be used for evaluation 
and assessment of the mathematical model. As time passed, a 
closed circuit cooling canal system was constructed for the 
Turkey Point facility and hot water was no longer discharged 
into the lower Biscayne Bay from this plant. The study plan 
was modified accordingly with the model study in Biscayne Bay 
concentrated on the Cutler Ridge Power Plant site. The first 
phase of the project covered the period March 1974 through 
October 1974. At the conclusion of the first phase study it 
was found that a division of the model into near and far 
fields would be more desirable to describe the details of the 
plume as well as the larger recipient domain. It was also 
found that in a large number of locations, namely inland 
lakes, surface waves were not an important factor and could 
be ignored to save computation time. This prompted division 
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of the modeling effort to rigid-lid and fr -\s-surface 
applications. The final report for the Phase I study 
includes the detailed rational and formulation of the 
separate programs. It should be noted that there are 
essentially two basic programs, the rigid-lid and the free- 
surface. They are applied to near and far field studies 
by appropriate modification to the boundary inputs to the 
model and the grid size to be used. However, the initial 
modification from far-field to near-field for either the 
rigid-lid or the free surface model is not trivial. In 
that sense there are four modeling efforts. During the 
period of the first phase study, the equipment for in-situ 
measurements was ordered, calibrated and tested. A data 
gathering test was carried out at the Biscayne Bay site. 

The theoretical study was initiated with a rigid-lid model 
(assuming fixed surface level). The model was formulated, 
the computer program developed and de-bugged. 

The present contract, which was initiated on October 15, 
1974 and covers a period of twelve months, covers the second 
phase of the model development. This study (refs. 3 and 4) 
serves the following purposes: (1) the development and 

testing of a generalized mathematical model to predict three- 
dimensional temperature and salinity distribution in coastal 
regions receiving hot discharges; and (2) the investigation 
of the thermal remote sensing systems by directly relating 
them to thermal radiation from the sea surfaces and by better 
accounting for the absorption in the atmosphere. This study 
program is being carried out in four parts. They are closely 
related and are to be conducted concurrently. The mathe- 
matical model development is the major part of the project. 
Remote sensing and in-situ measurements are needed to support 
the model development. At the same time, a turbidity remote 
sensing system is being developed. The relationship between 
these various phases of the stiidy may be shown by a flow 
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diagram in Fig. I-l. 

In this twelve-month period, the rigid-lid model has 
been developed and its application to far-field study has 
been completed. The velocity and temperature fields have 
been computed for different atmospheric conditions and for 
different boundary currents produced by tidal effects. The 
computations have been carried out for different time periods 
between 1 and 6 hours of real time. In parallel with the 
theoretical work, six experimental studies of the two sites 
in question (Biscayne Bay site and Hutchinson Island site) 
have been carried out. Four of the data gathering experiments 
were carried out at the Biscayne Bay site. Three of the ex- 
periments included airborne infrared data gathering by NASA-6 
aircraft. Two experiments have been carried out at the Hutch- 
inson Island site, one of which included airborne infrared 
data recording. The temperature fields obtained during the 
tests at the Biscayne Bay site have been compared with the 
predictions of the rigid-lid model study. The comparison of 
the theoretical and experimental results is encouraging. 

In parallel with the above described work, the development 
of a free-surface model has been initiated. In this model the 
water surface level changes will be included. The governing 
equations have been formulated and the computer programs have 
been written. De-bugging of the computer program is in progress. 

When the project is concluded, two important benefits will 
be obtained, viz., (1) it will be possible to select the proper 
sites for the future power stations so that they may cause 
little or no pollution, and (2) it will be possible to remotely 
monitor the discharges from the existing and future power plants 
to make sure that they do not exceed the allowable temperature 
limits and do not pollute . 
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II. GROUND TRUTH AND IN-SITU MEASUREMENTS 

Ground truth is primarily a service function for the other 
parts of the project, such as Model Development and Remote 
Sensing, and, consequently, is not an "end product" hy itself. 
Data from the field tests will be presented both here and in 
the Remote Sensing section. 

A. Background 

The major functions of this part of the project are 
to provide: 

1) Ground truth for evaluation of remotely 
sensed measurements, 

2) Boundary conditions and verification of 
the mathematical model, and 

3) Assistance in the development of new 
techniques for remote sensing. 

To perform these functions, the following duties are 

involved: 

1) Ordering of instrumentation, 

2) Calibration and testing of instruments, 

3) Assistance in coordinating and planning 

field tests (boat trips with and without aircraft) 

4) Supervision of measurements during field trips, 
and 

5) Reduction of data and preparation of data 
records . 

The quantities of interest in the measurement program 

are : 

1) Temperatures at the surface and at various 
depths measured by direct means , 
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2) Surface temperature measured by infrared 
sensor (heat spy) for comparison with 
airborne thermal scanner, 

3) Current speed and direction, 

4) Salinity, and 

5) Pertinent meteorological data, such as 
atmospheric temperature, wind speed and 
direction, humidity, and tidal information. 

Direct temperatures are sensed by thermistors and 
the accuracy of these measurements is limited only by the 
calibration. At present, this is about 0.2°C, which is more 
than adequate both for correlating with the remote readings 
and providing the necessary information for the mathematical 
model. The other measurements will be discussed elsewhere. 

B. Measurement System for Biscavne Bay 

Most of the instruments discussed in detail in previous 
reports are currently being used for this phase of the project. 
Several changes have been made, however, to correct certain prob- 
lem areas and to improve the system. Measurement of the re- 
sistance of the thermistors was somewhat difficult because the 
boat motion deflected the galvanometer in the Wheatstone bridge. 

A digital type null device was designed to correct this fault 
and it operated satisfactorily. However, the bridge-null indi- 
cator measurement is still a fairly slow procedure. A direct 
resistance indicator with digital display (Fluke Model 8000^) 
has been used during the recent tests and has performed well. 

The accuracy is consistent with system requirements and the 
measurements have been repeatable. 

As stated in previous reports, the modification to 
the General Oceanics Current Meter, which provides for audible 
signals to indicate the number of revolutions of the impeller 
blade, has proven very satisfactory for averaging the current 
flow over a short period, thus resulting in an average value 



_9_ 


rather than an instantaneous one which could be misleading. 

The counter also enables us to measure very low currents, which 
are prevalent in many parts of Biscayne Bay. 

The salinity measurements have been a problem. 

Direct conductivity type instruments have proven unsatisfactory 
due to the critical nature of the platinization of the elec- 
trodes. This appears to be a common fault of this type of in- 
strument’/ accordingly, an induction type meter was purchased 
(Beckman Type RSS-3) , which measures conductivity and hence 
salinity by noting the effect the conducting liquid has on a 
calibrated oscillator. This also has the advantage that it 
can be quickly recalibrated by inserting a known resistance in 
a loop through the sensor. The initial calibration was done 
chemically using known samples of saline solutions. Performance 
appeared to be satisfactory during the April 15 field trip but, 
unfortunately, several readings were beyond the scale of the 
instrument due to the highly saline condition of portions of 
Biscayne Bay, a result of the prolonged dry spell in this area. 
Satisfactory readings were obtained during the Sept. 4, 1975, 
field trip in Biscayne Bay and the Sept. 18, 1975, measurements 
off Hutchinson Island. 

G. Instruments for Hutchinson Island 

Several additional instruments were added for the 
experiments off Hutchinson Island because the water at the out- 
let of the nuclear power plant is much deeper. Rather than 
splice the wires of the present thermistors, several new tem- 
perature sensors were purchased with 60 feet of cable. Under- 
water splices may not be reliable and could be a serious problem 
if a failure occurred during a field trip. 

A new current meter was also purchased for these mea- 
surements with pressurization for a depth of at least 100 feet. 
The new instrument also has a built-in compass to give current 
direction, a pressure indicator to read depth, and a temperature 
sensor (thermistor) to give the temperature (within 1 C) where 
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the current is measured. This unit, the Endeco Type 110, 
has been tested and appears to operate satisfactorily. For 
low currents, the direct reading of the meter is somewhat 
inaccurate since the current flow is not steady. A counter 
on the number of revolutions of the impeller blade for a 
given unit of time allows the current to be averaged. 
Accordingly, a counter unit has been built with a self-con- 
tained timer with intervals of h, 1, and 2 minutes. The 
counter has a maximum rate of 1000 counts per minute, which 
can handle currents greater than four knots since the rotor 
unit develops pulses whose calibration is 3.56 pulses per 
second per knot . 

Other instruments for the measurement of resistance, 
salinity, atmospheric conditions, etc., will be the same as for 
the Biscayne Bay tests. 

D. Field Trips 

During the period of this contract, from Oct. 15, 

1974, to Oct. 15, 1975, tests were conducted as follows: 

Biscayne Bay: Oct. 24, 1974 

Apr. 1, 1975 
Apr. 15, 1975 
Sept. 4, 1975 

Hutchinson Island: May 1, 1975 

Sept . 18 , 1975 

During the previous contract year, tests were per- 
formed in Biscayne Bay on July 29, 1974, Aug. 21, 1974, and 
Sept. 24, 1974. The data from these tests were reported in the 
final report for the previous contract - 

The Oct. 24 test was conducted in the central region 
of Biscayne Bay in conjunction with the NASA 6 aircraft and the 
data sheets are included in the appendix. The major measure- 
ments on this data were the direct temperatures, the currents, 
and the surface temperature using an infrared sensor. A com- 
parison of the temperatures obtained from the ground truth and 
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the aircraft is given in the Remote Sensing Section, and the 
current information has been used in the Mathematical Model, 

The Apr. 1, Apr. 15, and Sept. 4 field trips 
focused on temperature measurements in the thermal plvime from 
the discharge canal of the Cutler Ridge Power Plant, a fossil 
fueled facility. Although this type of generating unit does 
not have as large an increase in water temperature as a nuclear 
power plant, it is nevertheless quite appreciable and very use- 
ful for testing the mathematical model. The morning measure- 
ments during the Apr. 15 trip were in the central region of 
Biscayne Bay and the data is compared with the aerial results 
in the Remote Sensing Section. A single measurement in the 
central region of the Bay was taken on Sept. 4 in the after- 
noon for calibration purposes. Data sheets for all these tests 
are given in the appendix. - 

Maps of the near field data location points for the 
Apr. 15 and Sept.! 4 trips are given in Figs. II-l and II-2. 

The time of measurement and surface temperatures are also 
provided to facilitate reference to the data sheets,. The sur- 
face temperatures are indicated on the isotherms obtained from 
the aerial scanner so that measurement comparisons can be more 
easily noted. These latter diagrams are given in the Remote 
Sensing Section. Isotherms for the tests on these dates are 
given in Figs. H-3 and II-4 for the extension of the center- 
line of the discharge canal into the Bay. Here the isotherms 
are shown where the depth and the distance from the mouth are 
the coordinates. This is useful in indicating the dissipation 
of the hot water as it is discharged into the Bay and could be 
compared with that predicted from the model. 

The Apr. 1 test was conducted primarily for famil- 
iarization with that region of the Bay, including navigation 
and techniques, for boat handling in very shallow water* Many 
surface measurements were made along the canal centerline 
although the navigational readings of some of the points are 
apparently in error. A diagram of the surface temperature as 
a function of distance from the mouth is presented in Fig. II-5 
and this could be of use in analyzing other data and examining 
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the results from the mathematical model. Similar diagrams 
for the Apr. 15 and Sept. 4 data are illustrated in Figs. 

H-6 and II-7. 

A preliminary test was made on May 1, 1975, off the 
coast of Hutchinson Island in the vicinity of the future site 
of the power plant outlet. To determine the effect of the 
outlet water, it would be advisable to measure conditions at 
several different times of the year before the plant goes into 
operation. The major purposes of this test were to find suit- 
able landmarks for navigation and location, examine temperatures 
at various depths, and record currents. 

For navigation purposes, the landmarks which would 
prove most useful are the power plant itself, a radio- television 
tower north of the power plant, and a small concrete tower 
south of the power plant. There is a dock near the outlet pipe 
which could be used for location. The inlet pipe area is 
marked by buoys and these would also be helpful. The system 
of navigation recommended here should include the depth. For 
example, the 18 ft. depth contour runs almost parallel to the 
coast and the boat's depth sounder can place us exactly at a 
given depth; other depths follow similar contours. The other 
dimension can be determined by sightings on the landm^j,'ks 
mentioned above. 

A test of conditions during the latter part of summer 
was conducted on Sept. 18, 1975, off the coast of Hutchinson 
Island. Unfortunately, weather conditions were poor and the boat 
trip had to be cancelled after the morning runs. Several data 
locations were obtained, however, and the data sheets are in- 
cluded in the appendix. Several interesting observations can 
be made from the data. First, the salinity was constant verti- 
cally and this is of value in deciding how much a factor 
salinity is in the model. The currents at all depths flowed 
approximately northward as is to be expected and this will be 
of considerable value in the model studies. The vertical 
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temperature gradients were slight, but this is normal for 
this area. All instruments functioned as expected and the 
techniques for taking ground truth measurements were improved 
as a result of the experience. 

The nuclear power plant will not be in fpll operation 
until early 1976, an,d present plans call for oceanographic 
measurements to be made in this vicinity during the fall and 
early winter of 1975. When the plant is in operation, simi- 
lar measurements will be acquired and the effect of the plant 
can be noted. More specific plans including test dates will 
be determined at a later time when availability of the aircraft 
and other factors can be taken into account. 
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III. REMOTE SENSING WATER SURFACE TEMPERATURE 

Four thermal IR remote sensing experiments have been 
conducted in Biscayne Bay iising NOAA-3 and NOAA-4 polar- 
orbit satellites and the Kennedy Space Center NASA-6 Beech- 
craft C-45H aircraft. These were conducted on July 29, 1974, 
October 24, 1974, April 15, 1975, and September 4, 1975. 

Isotherms depicting a warm water plume from Cutler Ridge power 
plant are presented for July 29, 1974 and April 15, 1975 from 
the 8-14 micron thermal IR channel of the Daedalus multispectral 
scanner aboard the NASA/KSC aircraft. Similar aircraft data 
are being processed for September 4, 1975. A weak isotherm 
pattern was mapped in Biscayne Bay on October 24, 1974 because 
the Cutler Ridge power plant was not operating and so did not 
produce a plume. 

Digital maps of remotely sensed temperatures for Biscayne 
Bay and surrounding areas from NOAA-3 and NOAA-4 satellites 
using 10,5 to 12.5-micron IR sensors are presented for October 
24, 1974 and April 15, 1975. Similar data are being processed 
for September 4, 1975. On July 29, 1974, extensive cloud cover 
prevented the NOAA-3 satellite -borne sensors from seeing Bis- 
cayne Bay at the time of the morning overpass. However, the 
aircraft was able to obtain data by flying beneath the cloud 
base. 

Unfortunately, the spatial resolution of thermal IR 
measuring systems used in the present generation of satellites 
does not provide as much detail as is desirable for this project, 
although the spatial resolution of these thermal sensors is 
improving with each new generation of satellites. Clouds may 
prevent satellite-borne sensors from seeing the ocean surface 
on some occasions but the frequency of observations is still 
better than that obtainable by other means for similar costs. 

It is believed that, ultimately, satellites can be used not only 
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for global ocean surface temperature mapping but should be able 
to provide the accurate, high-resolution water temperature data 
needed for monitoring thermal pollution from power plant dis- 
charges, detecting cold water upwelling and springs, and for 
model studies. 

A. Corrections to Satellite and Aircraft Measurements 
of Sea Surface Temperature 

Corrections must be applied to the satellite and air- 
craft measurements in order to obtain the real sea surface 
temperatures. These corrections are principally: 

1) departures of the sea from a true blackbody 
condition, and 

2) effects of the atmosphere through which the 
electromagnetic waves pass (Hiser et al.)Cl]. 

Investigations have shown that the correction for 
departure from a true blackbody is small [2]. It is probably 
within the limits of instrumentation error for our purposes 
(less than 0.5°C) and certainly much less than the atmospheric 
effects. A radiometer measures total radiance, R, which is the 
flux of radiation per unit solid angle per unit wavelength per 
unit area. The total radiance is composed of an emitted and a 
reflected component, namely: 

R = EB + rN (3-1) 

tJ'~re R is the total radiance of a surface, B is the radiance of 
a blackbody at the surface temperature, E and r are the emittance 
and reflectance of the surface, respectively, and N is the total 
radiance falling on the surface (the sky radiance). The re- 
flectance term in equation 3-1 is quite small for the ocean sur- 
face except when a scanning radiometer looks directly at a small 
sun igliht spot on a ^relatively smooth sea. Therefore, the 
radiance R measured by a radiometer would be nearly that of a 
blackbody at the sea temperature (EB) if there were no trans- 
mission losses through the atmosphere. 


ABSORPTION 


4 




The infrared region of the spectrum from 8 to 14 
microns, which is suitable for measuring terrestial tempera- 
tures of the order of 300°K, is subject to some absorption 
losses by gases and water vapor in the atmosphere (Fig. XII-1) . 



Fig. III-l. Absorption curves of H^O, CO 2 , and O^, 

between 3 and 14i|i . 

Ozone (0^) has a strong absorption band between 9 and 10 
microns; however, it is primarily in the upper atmosphere and 
does not affect data gathered by aircraft. The VHRR on the 
NOAA satellites operates in the 10.5 to 12 ; 5-micron region 
which is outside this ozone absorption band. 

Water vapor is the majoi variable in this thermal 
IR region of the spectrum with h minor contribution by CO^. 

If N is the radiance of the ocfean surface observed at zero 

..'S'- i' . 

altitude and N is the radiance as observed in a narrow spectral 
region from a satellite, the correction to the satellite ob- 
served radiance for the intervening air column is given by: 

-N= - J (Ns - Ijdu (3-2) 

o : 

Here B is the radiance of a blackbody at air temperature, t ( u) 

Si , 

is the transmittance of a mass of absorbing gas u measured 
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downward from the sensor, and U is the absorbing mass (mostly 
water vapor) for the entire column between sensor and surface. 

It is possible to calculate the amount of atmospheric 
absorption if the distribution of absorbing gases in the atmo- 
sphejre is known. The problem is complicated because the absorp- 
tion varies with pressure and temperature in a non-linear fash- 
ion. The calculations utilize laboratory measurements of 
absorption by atmospheric gases. However, because environmental 
conditions in the laboratory experiments are different from those 
in the atmosphere, theory must be applied. Computations of the 
outgoing radiation for 106 model atmospheres have been made by 
Wark et al. [3]. For clear sky conditions, water vapor was the 
primary variable. A correction for an average ozone distribu- 
tion was made; the error due to different ozone distributions 
is likely to be relatively small. 

The results for clear sky conditions and nadir angles 
of 0^ and 45° are shown in Fig. III-2 [4]. In this figure, the 
ordinate is the amount of precipitable water in the atmosphere. 
Horizontal lines indicate "average" conditions (W=2.25 cm) and 
"moist" conditions (W=3.9 cm). The nearly vertical lines are 
labeled T^, the equivalent blackbody radiation as measured by 
the satellite. Thus, for a 0° nadir angle with T^ = 290°, the 
correction to be added to obtain the temperature of the radiating 
surface is 6.8°C for average humidity and 9°C for moist con- 
ditions; at a 45° angle the corrections are 9.5° and 12°C, 
respectively. 

The precipitable water (W) for application in Fig. 

III-2 can be obtained from radiosonde data in many places. 
Radiosonde balloons are launched daily throughout the world at 
OG and 12 Greenwich Mean Time (GMT or Z). They obtain pressure, 
temperature, and relative humidity data to altitudes in excess 
of 50,000 ft. Most of the atmospheric moisture is contained in 
the troposphere below 30,000 ft. Therefore, these soundings 
can be used to calculate W for the atmospheric column. 


j 


5 6 

AT (®K) 


AT (*K) 


Fig. III-2. Clear-Sky Temperature Corrections 
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If the sounding is plotted on a Pseudo-Adiahatic 
chart, the mean vapor pressure (e) in millibars can be 
determined graphically for atmospheric layers between plotted 
points on the sounding. The water content for each of these 
layers can be calculated by using the equation: 

W cm = e In 

where g is the acceleration of gravity, is the pressure at 
the top of the layer in millibars and P 2 is the pressure at 
the bottom in millibars [s]. The total W is then obtained by 
summing these values for all layers of the sounding. Table 
III-l gives the precipitable water computed from the Miami air- 
port sounding at 1200Z (0800 EDT) on July 29, 1974. The right 
hand column gives the accumulated total W from the surface up- 
ward. The highest value of W in the top right corner is used 
to enter the diagram in Fig. III-2. Table III-l also provides 
data for lower altitudes that can be used to apply corrections 
to IR temperature measurements made from aircraft. 

B. Aircraft Data, July 29, 1974 

Fig. III-3 shows the surface isotherm pattern for a 
well developed thermal plume in Biscayne Bay from the Cutler 
Ridge Power Plant at 1230 EDT on July 29, 1974. These iso- 
therms were constructed from the calibrated gray-scale, 70-mm 
film transparency that was prepared from the magnetic tape used 
aboard the KSC NASA-6 aircraft to store the analog readout of 
the Daedalus multispectral scanner. Other data were taken at 
0930, 1030, and 1130 EDT, but the pattern at 1230 was most 
evident because the plume had time to develop as the power plant 
output increased and the discharge temperature reached a maximum 
The aircraft flew at 1500 ft. altitude and made twelve west to 
east data passes across Biscayne Bay between 0800 EDT (1200Z) 
and 1300 EDT (1700Z) . No ground truth data are available within 
the plume at the time the aircraft pass was flown. 
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TABLE III-l. PRECIPITABLE WATER AT MIAMI, 1200 Z (0800 EOT), 
July 29, 1974 


Approximate 
Altitude Feet 

Layers Pj^/P2 

Mean e for Layer 

Ciunulative 
W cm for Layers 

23,600 


400 

1.38 

5.376 

21.400 


438 





438 

2.25 

5.299 

18 . 300 


500 



■ 


500 

3.1 

5.113 

17,400 

* 

518 





518 

3.7 

5.045 

16.000 


550 





550 

4.4 

4.907 

14.600 


580 



' 


580 

5.1 

4.762 

13, 800 


598 





598 



■■ V 



6.2 

4.665 

11.400 


658 





658 

.,7.0 

4.296 

9,900 


700 





700 

9.8 

4.027 

7,100 _ 


780 





780 

10.6 

3.367 

6.700 


790 





790 






12.4 

3.283 

4,800 


850 





850 

17.8 

•2.718 

‘3,300 


895 





895 

22.8 

2.1479 

2.000 


940 





940 

29.0 

1.4516 

700 


985 





985 

28.0 

.6081 

. Surface 


1020 























Surface Tsothenns 

Cutler Ridge Power Plant , 

1230 EDT, Fit. Alt. 1500 Ft. i 
July 29, 1974 

Ground Truth Measurement, 2.5 NM (4,630m) 
South of Canal Mouth at 1200 EDT 
was 34.1°C = 93.4“F 

240 m. 


Cupola 

Fig. III-3 
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C. Satellite Data, October 24, 1974 

Temperature measurements of the earth's surface made 
by IR from satellites must be corrected for water vapor ab- 
sorption in the atmosphere. The Miami radiosonde data for 
0800 EDT on October 24 , 1974 has been used to determine the 
precipitable water in the atmosphere for this study in Bis- 
cayne Bay (see Table III-2) . The diagram in Fig. III-2 is 
used to determine the correction (At) to be added to the 
satellite measurements of sea surface temperature (T^,) in the 
presence of a given amount of precipitable water (W) . This 
diagram gives values for use near the satellite subpoint (Nadir 
Angle = 0°) and at 45° from the subpoint (Nadir Angle = 45°) . 
Linear interpolations for intermediate angles are within the 
limits of accuracy of the measurements. 

Figs.. III-4 and III-5 are digitized and contoured 
temperature maps prepared by computer from the NOAA-3 VHRR 
(10.5 to 12.5 microns) IR data in the vicinity of Miami on 
October 24, 1974. Fig. III-4 is a 1-km grid based upon the 
minimum resolution of the satellite and Fig. 3-5 is a 2-km 
grid with space -averaged data plotted at each grid point. The 
coldest areas (stippled) are radiation from clouds and can be 
ignored for our purposes. The satellite-measured isotherms in 
Fig. III-4 are 20°C (solid) and 18°C (Dashed) . In Fig. IXI-5 
they are 20°C (solid) and 24°C (dashed) . The direction of north 
is shown on each of these illustrations because the NOAA-3 
satellite track is inclined 16° from the meridians. 

The data in both of these figures are uncorrected 
for water vapor absorption. The Miami radiosonde at 1200Z 
(0800 EDT) on October 24, 1974, which was less than two hours 
before the satellite pass, gave a total liquid water content 
of 4.09 cm (see Table III-2) . Since the satellite measured 
surface temperatures of interest are near 290°K and the Nadir 
Angle is near 0°, the correction to be added is 9°C (see 
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Table III-2. Precipitable Water at Miami, 

1200 Z (0800 EDT), October 24, 1974. 


Approximate 
Altitude iFejet 

Layers 

Mean e 
for Layer 

Cumulative 
W cm for Layers 

23,600 

400 

0.16 

4.09 

18.300 

500 




500 

0.68 

4.05 

16.900 

530 




530 

0.7 

4.01 

9,900 

700 




700 

0.9 

3.82 

8.100 

750 




750 

3.1 

3.76 

7,800 

760 




760 

8.5 

3.71 

6.400 

795 




795 

9.8 

3.33 

4.800 

850 




850 

12.6 

2.68 

3,000 

910 




910 

15.8 

1.82 

200 

1000 




1000 

16.5 

0.33 

Surface 

1020 








































Fig. III-4. 


Digiti?.ed IR Temperatures from NOAA-3 
Satellite, October 24, 1974, 1342 Z 
(0942 EDT), 1-ltm Grid Squares. 
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Fig. III-2) In Figs. III-4 and III-5, the solid line is the 
20°C indicated (29°C water-vapor corrected) isotherm, and the 
dashed isotherm in Fig. III-5 is 24°c indicated or 33°C 
corrected. The temperature resolution of the VHRR sensor is 
approximately io.5*^C for a 300°K scene. This can cause the 
computer to produce some small contoured areas that may not be 
significant. The area of these is usually one resolution cell 
of the sensor and they are indicated to be only one degree 
colder or warmer than their surroundings. Most of these small 
contoured areas are not visible in the reproduction of Fig. 

III-4 and the computer averaging in Fig. III-5 eliminates them. 

These computer printouts were prepared by the 
National Environmental Satellite Service of NOAA in Washington, 

D. C. Because of the great enlargement of a small area for 
mesoscale analysis, a latitude and longitude grid coordinate 
system was not available in the program for these printouts. 

They must be located geographically by use of photographs made 
from the visual band of the multi-spectral scanner aboard the 
satellite. Such photographs were provided along with the digi- 
tal printouts for this purpose. A best-fit map of Biscayne Bay 
has been sketched on Fig. III-4, and a map of the southeast 
Florida coastal area with some latitude and longitude coordinates 
has been sketched on Fig. III-5. 

Both of these illustrations show slightly warmer 
water in Biscayne Bay than on the continental shelf, and Fig. 
III-5 shows warmer water in the Gulf Stream than on the shelf. 

The warmest water is shown in Fig. III-5 in the shallows north 
of Bimini and east of the Gulf Stream. The metropolitan area 
of Miami is shown slightly warmer than the glades and farm land 
to the south and southwest of the city. An interesting feature 
of these illustrations is the cool tongue of water extending 
into the central part of Biscayne Bay from the shelf area just 
north of Elliot Key. Low tide at Government Cut at the south 
end of Miami Beach occurred at 0950 EDT on this date. Thus, 
tidal currents should have been near a minimum at the time of 
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the satellite pass. More observations are needed to determine 
the significance of this feature. 

D. Aircraft and Satellite Data. April 15, 1975 

Fig. III-6 shows the thermal IR data gathered along 
flight paths 10 through 16 on the morning of April 15, 

1975 by the KSC NASA-6 aircraft. Clouds kept the aircraft at 
lower than the optimum altitudes for good aerial coverage by 
the thermal scanner, which has a 76-degree field of view nor- 
mal to the flight path. Boat measurements of ground truth 
temperature data were taken at points A, B, C, and D. The 
isotherms shown are uncorrected for atmospheric water vapor 
absorption. Since the theoretical corrections for water 
vapor are too large in some instances, the boat measurements 
are used to correct the aircraft data before final analysis. 
This figure illustrates that the aircraft thermal IR data 
are not synoptic in time and that the aerial coverage is in- 
complete because of cloud limitations on flight altitude. 

Fig. III-7 shows the thermal plume from Cutler Ridge 
Power Plant at 1546-1547 EDT as measured with the multispectral 
scanner (8-14 microns) from 2,000 ft. on flight line number 16. 
The highest indicated temperature in the plume is 89°F. The 
boat measurement at point A on the 83° isotherm, 31 minutes 
before flight time, indicates a +2.5°F correction to be applied 
to the IR data. This assumes steady state conditions in the 
plume for the 30 -minute period. Fig. III-2, using 2.07 cm 
of precipi table water for the flight altitude of 2,000 ft. from 
Table III-3, gives a much larger correction, which is probably 
excessive. Therefore, the +2.5°F ground truth correction was 
applied to the aircraft data and the adjusted isotherm values 
are shown in parentheses. 

An interesting feature of this plume can be noted by 
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Table III -3. Precipitable Water at Miami, 

1200 Z (0800 EDT), April 15, 1975 


Approximate 
Altitude Feet 

Layers ^ 2 

Mean e 
for Layer 

Cumulative 
W cm for Layers 

23,600 

400 

0.64 

5.15 

22,200 

425 




425 

0.95 

5.11 

20,300 

460 




460 

1.22 

5.03 

19,000 

485 




485 

0.56 

4.97 

16,600 

535 




535 

0.60 

4.91 

9,800 

700 




700 

1.22 

4.75 

9 , 300 

715 




715 

1.92 

4.73 

6,900 

785 




785 

5.6 

4.55 

6,200 

806 




806 

11.5 

4.40 

4,800 

850 



■ 

850 

. 

13.9 

3.79 

4,5G0 

860 




860 

16.8 

3.63 

2,600 

920 




920 

. 

19.5 

2.49 

2,000 

940 




940 

24.4 

2.07 

200 

1000 




1000 
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comparing it with Fig. III-3 for July 29, 1974, The plume on 
July 29 was elongated north-south while the one on April 15 
has its major axis oriented east-west. The difference in 
wind and tide conditions are. believed to be responsible for 
this. On July 29, 1974, low tide occurred at the Cutler Ridge 
Plant at approximately 1300 EDT while on April 15, 1975 high 
tide occurred there at approximately 1230 EDT. The outgoing 
tide prior to the image on July 29 probably contributed to 
expansion of the plume away from the shore, while the incoming 
tide on April 15 contributed to lateral spreading along the 
coast. Another major contributing factor may have been the 
surface wind conditions. Around noon on April 15, the winds 
were from the south at 13-16 knots blowing against the dis- 
charge plume. On July 29, 1974, they were much lighter and 
more variable, south-southeast at 6 knots at noon and west- 
southwest at 3 knots by 1300 EDT. Such light and variable 
winds could allow or help the plume to extend away from shore. 

Fig. III-8 is similar to Fig. III-6 except that it 
is for the afternoon aircraft flight on April 15, 1975. The 
isotherms shown are uncorrected for atmospheric water vapor 
absorption. The next step in the analysis was to apply a cor- 
rection based upon boat measurements of ground truth tempera- 
tures. These results were then extrapolated to construct an 
isotherm map (Fig. IV-40 in the Mathematical Model section 
of this report) . 

Fig. III-9 is a digitized and contoured temperature 
map prepared by computer from the NOAA-4 VHRR (10.5 to 12.5 
microns) IR data for south Florida on April 15, 1975 at 1335Z 
(0935 EDT). The data points are on a 2-km grid. Each data 
point is the average of six pixels or resolution cells of the 
satellite. The coldest areas (shaded) are radiation from clouds 
and can be ignored for our purposes. The satellite-measured 
isotherms in Fig. III-9 are 24°C, 28°C, and 32°C (solid) and 
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30°C (dashed) . 

This computer printout, which was prepared by the 
National Environmental Satellite Service of NOAA in 
Washington, D.C., contains a "U. S. Standard Atmosphere" 
water vapor correction for IR absorption. The standard at- 
mosphere correction is approximately 5°C less than that ob- 
tained for conditions at Miami on April 15 as given by Table 
III-3 and Fig. III-2. Thus, an additional 5°C should be added 
to the temperatures in the Miami area in Fig. III-9 to cor- 
rect them for the additional moisture of this climate. Apply- 
ing this correction gives some small areas with 35°C water 
temperatures in Biscayne Bay. These agree favorably with 
some of the highest temperature measurements made by boat on 
April 15. However, the spatial resolution and accuracy of 
the NOAA-4 (ITOS series) satellite makes these observations 
of limited significance for mesoscale thermal mapping, as 
has been discussed earlier in this section. 

E . Aircraft and Satellite Data, September 4, 1975 

Fig. III-IO is an IR temperature measurement of 
the Cutler Ridge Power Plant thermal plume made from the 
NASA-6 aircraft on September 4, 1975 at 1525 EDT, This was 
made from a flight altitude of 1800 ft. The temperatures 
are in degrees Fahrenheit and are not corrected for atmos- 
pheric water vapor absorption between flight level and the 
water surface. The power plant had been in operation for 
several hours and the plume was well developed by this time. 

Table III-4 gives the precipi table water as 
measured by the radio-sonde at Miami airport, 0800 EDT on 
September 4 > 1975. This information together with Fig. III-2 
and ground truth is used tc apply water vapor corrections to 
the NASA-6 aircraft and NOAA-4 satellite radiometer measure- 
ments of water temperatures in Biscayne Bay on September 4, 1975 
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Figure IlI-lO. Thermal Plume from Cutler Ridge Power Plant, 
1525 EDT, September A, 1975, Flight Line 16A 
Flight Altitude 1800 ft. 
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Table III-4. Precipitable Water at Miami, 

1200 Z (0800 EDT), September 4, 1975 


Approximate 
Altitude Feet 

Layers 

Mean e 
for Layer 

Ciittiulative 
W cm for Layers 

22,400 

420 

1.24 

7.577 

18.300 

500 




500 

0.42 

7.361 

16.800 

530 




540 

0.72 

7.337 

15.900 

550 




550 

2.9 

7.324 

14,200 

590 





590 


7.120 

9.800 

700 




700 

10.6 

6.146 

4.800 

850 




850 

15.8 

4.088 

3,600 

885 




885 

16.2 

3.451 

3.200 

900 




900 

19.0 

3.179 

2,800 

915 




915 

26.0 

2.865 

200 

1000 




1000 

31.2 
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Fig. III-ll shows the thermal IR data gathered 
along flight paths on the afternoon of September 4 , 1975. 
Clouds kept the aircraft at lower than optimum altitudes 
for good aerial coverage by the thermal scanner which has 
a 76-degree field of view normal to the flight path. The 
temperatures in Fig. III-ll are uncorrected for atmospheric 
water vapor absorption. This figure is similar to III-6 
and likewise illustrates that the aircraft thermal IR data 
are not synoptic in time and that the aerial coverage is 
incomplete because of cloud limitations on flight altitude. 

The flight strip data from Fig. III-ll were inter- 
polated and extrapolated to produce the complete isotherm 
map in Fig. III-12. Boat measurements of ground truth 
temperature data were taken in the warm water plume near 
the Cutler Ridge Power Plant. Since the theoretical 
corrections for water vapor are too large in some instances, 
the boat measurements were used to adjust the temperatures 
in Fig. III-ll to obtain those in Fig. III-12. 

The NOAA-4 satellite IR temperature data in 
digital format have not been received for September 4 , 1975. 
An SMS-1/GOES satellite visual image at 1330Z (0930 EDT) 
shows scattered clouds over the southeast coastal region 
of Florida at the approximate time of the NOAA-4 satellite 
pass over this area. Thus it is believed that some N07U^-4 
thermal IR water temperature data for Biscayne Bay will be 
available after processing by NOAA/NESS in Washington, D.C. 
Later in the day, the cloud cover became more extensive and 
would have prevented such satellite measurements. 

Additional analysis will be performed on the 
September 4, 1975 data during the next phase of the project. 
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IV. MATHEMATICAL MODEL 
A. Introduction 

The need for mathematical modelling in predicting 
and monitoring thermal pollution vy/as discussed in previous 
reports by Veziroglu et al. (1973, 1974) . Predictive studies 
of ecosystems can only be made by mathematical models. A 
priori knowledge of the effects of disturbances is essential 
for environmental impact studies. Thus, the mathematical 
model is a crucial tool in decisions involving power plant 
siting, land development, etc. 

-The University of Miami team undertook to develop a 
methodology using remote sensing and numerical modelling to 
study thermal pollution. The need for remotely sensed data in 
modelling has been discussed by Sengupta et al (1974). The re- 
mote sensing effort has been discussed in detail in a previous 
section. In this section the development and verification of 
the numerical modelling effort will be presented. 

The Biscayne Bay and Hutchinson Island areas in 
Florida were chosen as sites for initial development and cali- 
bration of the models. Fig. lV-1 shows the general layout for 
Biscayne Bay. The Bay is approximately 25 kilometers long and - 
8 kilometers wide. It has several basins connected by shallow 
limestone sills. The Bay is isolated from the Atlantic Ocean 
by a series of keys. The major ocean efflux occurs through the 
safety valve area between Key Biscayne and Soldier Key. This 
is about 5 kilometers long with shallow limestone sills 1 to 
3 kilometers wide. There are about 20 channels in this area 
with a maximum depth of 5 meters. The lower part of the Bay 
has a number of creeks open to the ocean. The maximum depth 
in the Bay is 4 meters. The Featherbed Banks effectively di- 
vide the Bay into two parts. There are two Florida Power & 

Light power plants located on the shoreline. The pirn t at 
Cutler Ridge has an open cycle cooling system discharging into 
the Bay. The Turkey Point plant uses a series of cooling canals 
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Fig. IV-2 shows >the general location of the Hutchin- 
son Island site. Here the open ocean boundary conditions are 
dominant. The water at this site is considerably deeper than 
the Biscayne Bay site. The effects of surface height dis- 
placements are more significant. A nuclear plant with two 
units of 850 MW each will discharge heated water at this loca- 
tion. 

The hydrodynamics and thermodynamics of an ecosystem 
are controlled by geometry, meteorological conditions and 
physical characteristics of the water, such as salinity and 
turbidity. The governing equations are therefore composed of 
the Navier-S tokes equations, energy equation, salinity equa- 
tion, the equation of state and some transport equation for 
turbidity. Various assumptions can be made for different 
situations leading to simplification or elimin^.tion of equa- 
tions. 

The effects of salinity have been ignored in the 
models until now although capability to include them exists. 

A discussion on the equation of state in Appendix I clarifies 
the relative importance and role of salinity. The individual 
formulations discussed later in this chapter justify the in- 
clusion or exclusion of salinity more specifically. In general, 
if the salinity gradients are small enough to give densimetric 
Froude numbers greater than one, salinity effects can be ignored. 
Even where salinity does modify circulation, if the convective 
effects on temperature fields are non-significant, salinity can 
be ignored since our main concern is thermal pollution. 

In the initial feasibility report it was proposed to 
develop a single model to incorporate the effects of thermal 
pollution for the whole Bay. However, since the effect of the 
plume is manifest only in a small domain of the order of a few 
square miles, it is wasteful in terms of computational time to 
have a single model. Therefore, separate models are being 
developed for near field and far field. The near field is the 
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region affected by the plume. The far field is the whole eco- 
system under consideration. The far field affects the near 
field significantly but not vice versa. The initially sug- 
gested numerical scheme was a marker and cell method. In our 
present study it was felt that a stretching suggested by 
Phillips and used successfully in studies by Freeman et al.(1972) 
and Sengupta and Lick (1974) was more suitable. The method 
allows easy adaptation to various bottom topographies, an im- 
portant requirement for any general model. 

Owing to the aforementioned considerations, the 
mathematical modelling effort has been divided into two parts; 
the rigid-lid formulation and the free surface formulation. 

The rigid-lid formulation artificially forces vertical velocities 
to zero thereby eliminating surface gravity waves. This usually 
leads to larger integration time steps since the Courant-Levy- 
Fredrichs stability condition is relaxed. The effect of sur- 
face height variation manifests itself as a lid-pressure. The 
free-surface formulation allows surface movement and gravity 
waves, the surface pressure being atmospheric. The rigid-lid 
formulation is useful for lakes or cases where surface eleva- 
tion changes are slow and are small compared to total depth. 

The transients are somewhat distorted owing to damping. The 
free surface model is more general although usually its full 
capacities are utilized in situations where surface height 
fluctuations are important. 

No rigorous criterion exists for selecting either 
the free-surface or rigid-lid model. However, there are some 
comparative studies to guide the user’s judgment. Crowley 
(1969, 1970) compared both free-surface and rigid-lid models 
for ocean circulation and found the rigid lid model satisfactory 
except for errors in predicting Rossby wave speeds. The rigid- 
lid results for Biscayne Bay, when compared with free-surface 
results, will indicate the differences, if any, for shallow 
bays with tidal flushing. A conservative approach would be 
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to use the free surface model for all open ocean and tide 
dominated flows. Even then a vast majority of thermal 
pollution sites would be accurately modelled by the rigid- 
lid formulation. 

The rigid-lid far field model has been developed 
and applied to Biscayne Bay and general agreement with 
remote sensing data has been obtained. The free-surface 
model has been programmed and preliminary results obtained. 

B. Free Surface Model 

1) Model Formulation 

The set of equations governing the behavior 
of the fluid flow are those expressing the conservation of 
mass, momentum, and energy in a turbulent flow ^ The Cartesian 
coordinate system is used here, with z-coordinate in the verti- 
cal direction as shown in Fig. IV-3. In order to keep the 
generalized nature of the model, all the significant terms in 
the pertinent conservation equations are retained. Included 
are the effects of buoyancy, inertia, coriolis, density and 
turbulent mixing. Wind shear and heat flux at the free sur- 
face are also considered. Provision is made to include sal- 
inity effect by adding a salinity concentration equation and 
modifying the equation of state. The equation of state given 
in this section is in a general form. As discussed in the 
last section, a specified equation of state may be chosen 
depending on the local situation of the site where the model 
is to be applied. 

One major difficulty in the treatment of the 
Free Surface Model is at the free surface boundary. The 
boundary conditions can be specified but the position of the 
free surface is a variable and must be determined. In addition 
to this, the shape of the free surface is irregular, making 
it very difficult to apply any grid system at this boundary 
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for numerical solution. The approach used in this study 
is to follow a vertical stretching transformation suggested 
by Philips (1957). Using this transformation, the free surface 
becomes a fixed flat surface and the variable surface height 
remains constant. 

Another scheme used in this study is to adapt 
a horizontal stretching of coordinates similar to that of 
Waldrop and Farmer (1974). In this way, a variable size 
grid system may be used to give more detailed description 
of the flow in the thermal discharge area, while larger grid 
sizes may be used for farther points to cover a large area. 

It is obvious that closed-form analytical 
solution of this set of complicated, coupled, unsteady, 
three-dimensional non-linear partial differential equations 
is impossible to get. The finite difference method is used 
to obtain numerical solutions. The central difference 
approximation is used for the space differential. For the 
time differential, forward difference and central difference 
are used alternately to insure stability. The formulation 
of the free surface model is completed and has been pro- 
grammed for the UNIVAC 1106 computer system at the University 
of Miami Computer Center. The computer program is being 
debugged and some preliminary results have been obtained. 
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2) Governing Equations 

The following partial differential equations, written in 
Cartesian coordinates, describe the three-dimensional unsteady 
fluid flow: 

CONTINUITY EQUATION: 
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Here 

u, V, w are scalar components of the velocity vector 
in the X, y, z directions# respectively; 
p is the pressure; 
p is the density; 
f is the coriolis parameter; 
g is the gravity 
T is the temperature; 

S is the salinity concentration; 

Kjj and are the eddy viscosity in the horizontal and 
vertical directions; 

Bjj and are the eddy thermal diffusivity in the 
horizontal and vertical directions; and 
D„ and D are the eddy diffusion coefficient in the 

xl V 

horizontal and vertical directions. 

For the type of flows under consideration, transient response 
is gradual and the inertia and viscous terms in Equation (4-4) 
are negligible when compared with the pressure and the gravity 
terms. Equation (4-4) may then be simplified to give the 
hydrostatic approximation for pressure by 


0 = 


az ^ 


(4-8) 


3) Transformation of the Vertical Coordinate 

In the present study, both the free surface and 
the realistic bottom topography are considered, and to apply 
finite difference equations at these irregular boundaries for 
numerical solutions is quite difficult. To solve this problem, 
a transformation of the vertical coordinate is used to give a 
constant height for the free surface and a flat bottom boundary 
in the new coordinate system by letting 
a = X 


P = y 

Cf = Z {x,Yj z ,t ± ^ 

H(x,y, t) H(x,y, t) 


(4-9) 


where Z = Z(x,y,z,t) is the position of the fluid 
element relative to the free surface and 
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H = H(x,y,t) is the depth contour also relative 
to the free surface. 

The Z variable decomposes into the two components Z = z + 
'Tl(x,y,t), where z represents the vertical position relative 
to the mean water level (MWL) and T] gives the free surface 
elevation measured positively upward from this mean water level. 
The H variable separates into the components H = h(x,y) + 
Tl(x,y,t), where h represents the depth relative to the mean 
water level. The value of (T ranges monptonically from zero 
at the free surface to unity at the bottom of the lake. The 
relationship between these two coordinate systems is shown in 
Fig. IV-4. Fig. IV-5 gives the transformed vertical topography 
in the (a^P^cr) coordinate system. 

Substituting transformations (4-9) into Equations (4-1) to 
(4-8), the governing equations in the (a,P,a) coordinates may 
be obtained as follows: 

Continuity Equation: 

.an 


M + 5 ( .H u ) llHyl , 

at t ao + iJ + 


= 0 


Momentum Equations 

a (Hu) a (Huu) a (huv) ...a (uO) 
at ^ ba ap ^ ^ ha 

- i (|?) fv ] 

+ al- («!?) . ^ 


(4-10) 


P 

t Jcr (PKvsS) 3 

a(Hv) 

a (Huv) 

a (Hvv) 

at 

aa 

ap 


aa 


(4-11) 


= H c - i (|f) + g(g|f - fu ] 
‘Si 'sl ^ + 5! ^ 




(4-12) 


-56- 


Energy Equation : 
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Salinity Equation 
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The hydrostatic approximation for pressure: 

P(a) = P(0) + gH p (a) dcr 

Jo 

Here 0 = is the vertical velocity in the (a,p,a) coordinates 
It should be noted that n is zero at both the free surface and 
the bottom boundary. 

4) Horizontal stretching of the Coordinate System 

It is desirable to obtain a more detailed descrip- 
tion of the flow near the discharge point while larger grid size 
may be used in the further points to save computation time, A 
horizontal stretching of the coordinate system is used here to 
create a more efficient use of the grid points by letting 
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Here C^# and are arbitrary constants. Using this 

transformation, even increments of AX and Ay will produce a 
closer spacing of grid points in the original physical coordinates 
near the plume discharge where better resolution is needed. This 
transformation will also produce ever-increasing increments of 
the physical grid size away from the plume discharge, thus ex- 
tending the region of computation as shown in Figs. IV-6 and IV-7 . 

The transformed governing equations in the stretching (X,Y) 
coordinates are then r*iven below: 

Continuity Equation: 
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V - Momentum. Equation 
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The hydrostatic pressure approximation will remain the same: 


P(ff) = P(0) + gH J* p (a) da 
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5) Initial and Boundary Conditions 

To start the calculation, initial conditions 
are needed. Initial velocity, temperature, density and 
salinity concentration should be specified throughout the 
domain. This information may be obtained from remote sen- 
sing, in-si tu measurement and other sources. Boundary con- 
ditions are to be specified at the lateral boundaries of the 
domain, at the air-water interface and at the bottom of the 
basin. Also needed are boundary conditions at the efflux 
points and at the thermal discharges. At all solid boundaries, 
no-slip and non-normal velocity conditions are applied. 
Adiabatic conditions are assumed at these boundaries for the 
energy equation. At the air-water interface, the pressure is 
set equal to the atmospheric pressure. Also, wind stress and 
a heat transfer coefficient are to be specified at this boun- 
dary. At points of thermal discharge, both velocity and 
temperature are to be specified. 

6) Numerical Scheme and Solution Technique 
The finite difference equations for the set 

of governing equations are obtained by using finite difference 
approximations for the differentials. The central difference 
is used for space while the forward difference and central 
difference are used alternately for time to insure stability. 
The grid point system used in the calculation is shown in Fig. 
IV-9. A flow chart is given in Fig. IV-10 to show the pro- 
cedures used in the numerical calculation. The computer 
program is run on the UNIVAC 1106 Computer System at the 
University of Miami Computer Center. It is arranged to have 
each part of the calculation carried out by a subroutine. All 
of these subroutines are called in by the main program. To 
apply this model to a specific site, only those subroutines 
that involve the boundary and initial conditions need to be 
modified. All other subroutines are quite general. The 
programming for the free surface has been completed. 
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C. Rigid Lid Model 


4 


One of the major model development tasks proposed 
for the project was the development of a rigid- lid model that 
can predict the hydrodynamics and thermodynamics of basin wide 
phenomena. Specifically, the capabilities of such a model are: 

(i) Predict the wind driven circulation, 

(ii) Predict the circulation caused by inflows and 
outflows to the domain, 

(iii) Predict the thermal effects in the domain, and 

(iv) Combine the aforementioned processes. 

To this end a rigid-lid model has been developed and applied 
to Biscayne Bay. The calibration procedure consists of com- 
paring ground-truth corrected airborne radiometer data with 
surface isotherms predicted by the mathematical model. This 
procedure is complete and the model is ready to be tested at 
other sites. Lake Be lews in North Caroline, which Duke Power 
Company uses as a cooling lake, is a prospective site. 

Some of the results presented here have been previously 
discussed in contract reports by Veziroglu etal. (1974, 1975) 
to NASA-KSC. This section attempts to present a comprehensive 
summary of the application of the rigid-lid model to Biscayne 
Bay. This subsection will deal with the formulation and 
justification of a rigid- lid model, its solution, the results, 
and their verification with NASA-6 thermal scanner data. 

1) Rigid-Lid Formulation 

The governing equations for fluid flow and heat 
transfer in a domain are a system of coupled non-linear second- 
order three-dimensional partial differential equations, which 
satisfy the local conservation laws for mass, momentum and 
energy. An equation of state completes the set of equations. 

For laminar flows molecular transport properties can be used 
to describe heat and momentum transfer. However, most geo- 
physical flows are turbulent, so the time averaged Navier-Stokes 
equations and the turbulent closure condition are used. 
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in which the Reynolds stress terms are replaced by eddy 
transport coefficients. The governing equations under 
these assumptions for an incompressible fluid have been 
derived by Sengupta and Lick (1974) and have been presented 
in a report by Veziroglu et al. (1974) . The hydrostatic 
assumption has been made. The fluid is considered to be an 
incompressible Boussinesq fluid. These assumptions are well 
accepted and invoking them does not require further explanation. 

The significant assumption made here is the 
rigid- lid assumption. This assumption effectively eliminates 
surface-gravity waves by imposing a zero vertical velocity 
at the surface. This essentially means that the surface allows 
slip conditions without any normal velocity. The surface 
pressure is no longer atmospheric but represents a "lid- 
pressure", which under some special steady state conditions 
can be related to the surface elevation that would occur if 
no lid were present. Transient time scales are distorted by 
this assumption but circulation patterns are not affected 
greatly. This has been demonstrated by Berdahl (1970) , 

Crowley (1969, 1970), Hag and Lick (1973) and other investi- 
gators, For many environmental investigations, including the 
study of thermal pollution, the rigid-lid assumption is 
adequate. Considerable savings in computing time result from 
eliminating the necessity to reproduce surface gravity waves. 

The assumption is particularly applicable for cooling lake 
studies, where surface elevation changes are not large and do 
not occur rapidly. 

The salinity equation has been dropped since 
its effect in lakes is negligible and the rigid- lid model is 
particularly designed for lakes. Even for Biscayne Bay the 
effect of salinity on temperature fields is negligible, as is 
shown in the Appendix on the equation of state. However, 
salinity can be included in the model without too much difficulty. 

For computational simplicity vertical stretching 
is applied to the domain to convert it to constant depth basin. 

This has been used by Freeman et al. (1972) and modified for 
rigid- lid applications by Sengupta and Lick (1974) . The stretch- 
ing makes the model easily applicable to various bottom topographies i 
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Governinq Equations 

The programming difficulties for a three dimensional basin 
suggest a stretching of the vertical coordinate of the form. 


V = 


Vr (X,y) 


This converts the basin to a constant depth one. Now the same 
number of grid points in the vertical direction may be used at the 
shallow or deep parts of the basin, without using variable grid sizes. 
The details of derivation are presented by Sengupta and Lick (1974) . 
The non-dimensional governing equations in the new coordinate system 
ck', p, Y (instead of 5?',y,z)are: 
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Equation of State 

r = 1. 029431 -.000020T-. 0000048? 
where T is in degree celsius 


(4-29) 
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To obtain a predictive equation for pressure the horizontal 
momentum equations are integrated from z=0 to z=h where h is the 
non-dimensional depth h/fj . The integrated equations are then 
differentiated with respect to a and p and summed. 

The Poisson equation for surface pressure becomes: 
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The last term is the Hirt and Harlow (1964) correction term 

which accounts for non-zero vertical velocities at the rigid lid. 

The variables (B , B , and A . A etc.) are given below. 
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The set of equations (24-30) together with the appropriate 
boundary conditions constitute the mathematical model. 

Initial and Boundary Conditions 

The nature of the equations require initial and boundary 
conditions to be specif ied. As the initial condition the velocities, 
temperature and density are specified throughout the domain. 
Boundary conditions are specified at the air-water interface, geo- 
graphical boundaries of the domain, the bottom of the basin and 
at efflux points. At the air water interface a wind-stress and 
a heat transfer coefficient are specified. The conditions on the 
lateral walls are no slip and no normal velocity for the momentum 
equations. These walls are assumed to be adiabatic. At the 
floor of the basin, the conditions of no slip and no-normal 
velocity are applicable. The energy equation has a heat flux 
boundary condition, considered adiabatic for the present study. 

At points of efflux the velocities are specified and gradient 
of temperature normal to domain boundary is considered zero. 

These open boundary conditions are mos!: difficult to specify. 

Therefore, the boundary conditions are in summary:. 
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where and T^y are wind stresses in the x and y 
respectively 

At the bottom of the basin: 
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Method of Solution 

The solution technique involves explicit numerical solution 
of the governing equations. The predictive equations which 
predict velocities are the momentum equations. The energy 
equation is explicitly solved for temperature. The Poisson 
equation, derived from the vertically integrated horizontal 
momentum equations for surface or lid pressure, is the pre- 
dictive equation for pressure. 

The finite difference equations for momentvim and energy 
are obtained by using a forward time central space scheme with 
Dufort-Frankel differencing on the diffusion terms. 

The finite difference equation for the u momentum equation 
may be '-written as: 

u^^^-u^ = (Inertia)^ + (Pressure)^ + (Viscous)^' n+1, n 1 

At 

The convective stability criterion is still effective with 
some relaxation of the diffusive stability criterion. The 
pressure equation is approximated by a five point scheme and 
solved by successive over relaxation (Liebmann method) . 

A staggered mesh is produced in which horizontal velocities 
and temperatures are calculated at integral nodes and vertical 
velocity and pressure are calculated at half grid points. 

Solutions have been obtained for Biscayne Bay in three 
stages. First, a closed basin approximation was made and wind 
driven circulation was obtained. In the second step circulation 
in the Bay was obtained with an ocean efflux specified; the 
simultaneous effects of wind and ocean efflux were also investi- 
gated. The third and final stage obtained the temperature field 
for various ambient conditions and parameters. The results were 
compared with airborne thermal scanner data to calibrate the 
model. Some of the results have already been discussed in re- 
ports by Veziroglu et al. (1974). This procedure was followed 
in order to check the performance of the model for as wide a 
range of environmental situations as possible. 
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2) Wind Driven Circulation 

Wind driven circulations for a sustained wind 
of about 5 mph for various directions (with a corresponding 
wind stress of .1 dynes/cm ) have been obtained. Two cases* in 
which the wind is from the southeast and the north* are 
discussed in detail. 

Fig. IV- 10 shows the surface velocities for 
winds at 5 raph. The velocities vary significantly in the 
basin. The variation of velocity direction and magnitude 
results from a number of factors. The outline of the Bay 
guides the current near the shoreline. The bottom slopes 
tend to influence the velocity* and the local depth of the 
Bay is also significant. The shallower regions like shoals 
and near shore areas have smaller velocities. The reduction 
of velocity near the shore is caused mainly by the smaller 
depth rather than by shore boundary layers. Estimates by 
Sengupta and Lick (1974) have indicated that the sidewall 
boundary layers are thin for similar situations and do not 
extend as far as the nearest interior node (0.8 kilometers 
from the shore) . Fig. IV-11 shows the velocities at a depth 
of two meters below the surface. There are no velocity 
vectors in large parts of the domain since the Bay is less 
than two meters deep in most places. The return flow can be 
clearly seen. 

A cross-sectional view of the velocity field 
indicates a cell- like structure. Fig. IV-13 shows a section 
through J=7 (refer to Fig. IV-12) . It can be seen that the 
surface velocities are in the direction of the wind stress 
and sharp upwelling occurs in the near shore region. Wherever 
there is a sharp change in bottom depth* horizontal gradients 
of velocity become significant and require vertical velocities 
to satisfy the continuity equations. 

Fig. IV-13 shows the Featherbed regions (at 
point A) which divides ^the bay. In this area upwelling and 
downwelling regions are in close proximity. Due to this 
barrier (shallows) interchange between the different basins 
of the Bay may be greatly retarded. . 
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Wind from Southeast 

The prevailing winds during spring and summer are from 
the southeast. The surface velocities for a 5 mph wind from 
the southeast are shown in Fig. IV- 14. Most of the surface 
currents are in the direction of the wind. It is important 
to note the effect of the shoreline and bottom topography on 
the flow. The magnitude of the velocities varies with depth 
and is smaller near shore and in the shallow shoal regions. 

In some parts near the mainland, the velocities are almost 
parallel to the shore. Near Cutler Ridge (1=11, J=4) the 
velocities are towards northwest indicating that with a south- 
east wind the plume will be pushed back towards the discharge 
mouth; thus shortening its longitudinal extent. The return 
flow at a depth of 2 meters is seen in Fig. IV-15. A cell 
like flow pattern is evident in Fig. IV-16. The Featherbed 
Banks tend to divide the two regions quite effectively. 

3) Tidal Currents 

The rigid- lid model does not reproduce level 
variations due to the tidal cycle. However, the tidal driv- 
ing mechanism can be modelled by specifying the inflow and 
outflow velocities at ocean inlets during various phases of 
the tidal cycle. 

The current measurements available for Biscayne 
Bay are very sparse and not systematic. Some measurements for 
Angelfish Creek, Broad Creek and Caesar Creek do exist. From 
these measurements and an average estimate for the shoal area, 
an order of magnitude value for tidal efflux velocity has been 
obtained . 

A 10 cm/sec velocity for both the incoming and 
outgoing tidal phases is assumed. The program is executed 
with an incoming tide, and then the ocean-bay velocity pattern 
is gradually reversed to obtain a 10 cm/sec outgoing velocity. 
Though the details of phase relations, level changes and time 
dependence cannot be precisely modelled, the results should give 
a meaningful qualitative picture of the circuit, tion . 

Fig. IV-17 shows the surface velocities with an 
incoming tide. The major portion of the tidal mass influx 
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travels into the South Bay while the flow towards Ricken- 
hacker Causeway is minimal. The currents near the shoals 
follow the influx direction but turn southward near the 
shore. The water flow through the creeks causes local 
modification of this southward flow. The velocities are 
reduced around Featherbed Bank, causing relatively high 
velocities in the deeper regions around the Bank. The 
velocities also increase to the south as the Bay narrows. 

The current near Cutler Ridge is from east to west indica- 
ting that the plume (which issues north south) will be ■ 
turned towards the west shore of Biscayne Bay. Fig. IV- 18 
shows the velocities at a depth of 1 meter with an incoming 
tide. The flow is unidirectional at most points. The re- 
verse current seen in wind-driven circulation does not occur. 
The reason for this is that in the case of wind driven current, 
the driving mechanism is at the surface, resulting in a pres- 
sure gradient that causes reverse flow; whereas, in the case 
of tidal currents the driving mechanism is a flow at all 
depths in the ocean inlets. 

The incoming flow was gradually reversed. 

Fig. IV-19 shows the surface velocities at an intermediate 
state. The currents have reversed in some places but not in 
others. In Fig. IV-20, the currents at a depth of two meters 
are almost completely reversed. This phase lag between sur- 
face currents to currents at lower depths is expected; the 
surface currents are larger and so take longer to reverse. 

With outgoing tide from Biscayne Bay into the 
ocean the flows are as expected. Fig. IV-21 shows that the 
bulk of the outflow comes from south bay. At Cutler Ridge 
the current is now from west to east, therefore with outgoing 
tide the plume is expected to turn towards the east. Fig. 
IV-2la shows that the currents at a depth of one meter are 
completely reversed compared to the inflowing currents. 
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4) Currents with Wind and Tide 

To investigate the circulation pattern with 
wind and tide acting simultaneously, a combination of 10 
mph southeast and an incoming tide was considered. Fig. IV- 
22 shows the corresponding surface velocities for comparison 
with Fig. IV-17 (incoming tidal current only) . The effect 
of the southeast, wind is to turn the current vectors toward 
the northwest. This results in more water flow towards the 
northern part of the Bay. The velocities in the south bay 
are still essentially southward but are of a smaller magnitude 
due to the southeast wind. Figs. IV-23 and IV-24 show similar 
circulations at 1 and 2 meters depth, respectively. From 
Fig. IV-22 it can be seen that the current near Cutler Ridge 
is in a north -west-west direction indicating that the plume 
will be pushed toward the shore and bro:»dened. 

5) Temperature Fields 

The energy equation is coupled to the momentum 
equations through the bouyancy term. Temperature gradients 
in the domain cause density gradients and consequently density 
currents. The circulation modification due to density currents 
causes changes in the temperature field. It is of interest 
to investigate the role of density currents in modifying tempera 
ture distribution. This is also important from a computational 
point of view since weak coupling would justify solying the 
energy equation with an invariant velocity field calculated 
a priori. 

The temperature solutions for various combina- 
tions of parameters have been obtained for comparison to the 

t 

remote sensing data from the April 14, 1974 NASA-6 flight over 
Biscayne Bay. It should be noted that the remote sensing data 
is not synoptic. Table IV-I shows a list of some of the cases 
for which solutions were obtained. At first a number of cases 
were run to obtain an estimate of the surface heat, transfer 
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coefficient, the effect of wind direction, wind magnitude, 
tidal flow and air temperature* The Bay was assumed to be 
initially isothermal and results were obtained for heating 
during the mid morning hours. 

Cases 1 and 2 were run to investigate the 
effect of the surface heat transfer coefficient. The values 
for heat transfer coefficients are taken from Edinger and 
Geyer (1971) and suitably increased for daytime heating only. 
(The Edinger and Geyer values are averaged over the entire 
day.) A wind of 10 mph from the southeast is also assumed. 
Figs. r^^~25 and IV-26 show the surface isotherms after 3 hours 
of heating for Cases 1 and 2, respectively. The shallow near- 
shore regions are warmer. The Card Sound area appears as a 
somewhat separate basin with closed isotherms. Featherbed 
Bank shows up as a developing warmer region. The warm spot 
over these banks is more pronounced for Case 1 which has a 
higher heat transfer coefficient. 

Fig. IV-27 shows surface isotherms for Case 3, 
which are almost identical to those of Case 2. Increased 
convective effects due to higher wind shears has had little 
effect. 

Cases 4 and 5 were executed to investigate the 
effects of density variation. Figs.IV-28 and IV-29 show sur- 
face and sectional (vertical) isotherms for Case 4, where 
density effects are excluded. Fig. IV-30 shows the surface 
isotherms after 3 hours for the same case. Case 5 is identi- 
cal to Case 4 except it includes the effects of density. Sur- 
face and vertical isotherms for Case 5 are shown in Fig. IV- 
31 and Fig. IV-32. Comparing Figs. IV-31 with IV-28 and 
IV-32 with IV-29, the effects of density can be seen. Since 
warmer fluid tends to rise and spread, the warm spot oyer 
Featherbed Banks is larger for Case 5 than for Case 4. Also 
comparing vertical sections, some of the sharp isotherm 
gradients near boundaries have been reduced by density currents 
i.e., the isotherms have been made more horizontal. 
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Cases 6 through 10 were executed to investigate 
the effects of tidal flow. A combination of surface and 
initial temperatures and heat transfer coefficients were used 
to match the results with surface isotherms from remote 
sensing data. Figs. IV-33 through IV-39 show surface iso- 
therms and vertical sections for some of these cases. Case 8 
agreed most closely with the remote sensing data. Fig. IV-40, 
the surface isotherms obtained from IR data corrected with 
ground truth data, shows extremely good agreement with Fig. 
IV-35, the surface isotherms for Case 8. It should be pointed 
out, however, that the IR data is not truly synoptic. There 
is a time lag of almost three hours between the near shore 
flights and the flight over the keys. Considering this factor 
and the assumption of isothermal conditions in the early 
morning, the agreement is satisfactory. 

By considering the cases mentioned above a 
choice of parameters which might agree with IR data can be 
made. The measured quantities cannot be altered? namely, 
wind speed (10 mph southeast), incoming tide and ambient 
temperature (30'^C) . To minimize the effect of initial tem- 
perature conditions, the run should be executed with as early 
a starting time as possible. Case 11 shows the selected 
parameters. 

At 9 AM near location (1=11, J=3) the mea- 
sured temperature was 25.6°C with no vertical stratification. 
Since this is a near shore location the "average" temperature 
in the Bay can be assumed to be lower. Since the detailed 
temperature field was not known as an initial condition it 
was assumed that that Bay is isothermal at a temperature of 
24.5°C at 8 AM. This approximation in initial condition be- 
comes less important with time. Therefore, comparisons with 
ground-truth or IR data were made only after the program had 
been executed for a few hours. The comparison with mid- 
afternoon IR data was after 6 hours. 
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case 12 uses the same set of parameters hut 
includes effects of density. Figs. IV-41 and IV-42 show 
isotherms for Case 11 after 3 hours. Fig. IV-43 shows the 
surface isotherms after 6 hours, which may be compared with 
Fig. IV-40. The qualitative and quantitative agreement with 
the remotely sensed data is quite remarkable considering its 
non-synoptic nature and the difficulty in specifying the 
initial conditions for the Bay. Fig. IV-44 shows a transect 
for Case 11 after 6 hours. Figs. IV-45 and IV-46 show the 
isotherms after 3 hours with density effects included. 

Since they are only slightly different from Figs. IV-41 and 
IV-42, density currents are not significant in changing the 
temperature distribution during this period. 

The temperature distribution after 6 hours for 
Case 12 is shown in Figs. lV-47 and IV-48. It is evident 
from the surface isotherms that the hot spot over the shoals 
has increased in area. This is expected since warmer fluid 
rises and spreads at the surface when density effects are 
included. The colder region in the Card Sound area has de- 
creased in size as warmer fluid from the shallower shores 
has spread at the surface. The agreement with remote sensing 
data is very good, as a comparison with Fig. IV-40 indicates. 

The ground truth data were obtained only in 
the morning on April 15, 1975 field experiment. 

Results for the model for Case 12 after 3 hours 
of heating were compared with ground truth measurements at 
location (1=11, J=7) at 11 AM. Fig. IV-49 shows a transect 
of the bay along J=7 . The calculated isotherms are shown with 
ground truth data at location marked by asterisks. It can be 
seen that the agreement is to within 0.7^C. This is satis- 
factory for most environmental models. 
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6) Conclusions 

A three-dimensional rigid- lid model has been 
used to study the general circulation and natural temperature 
distribution in Biscayne Bay. Remote sensing data has been 
used to calibrate and verify the model. The relevant con- 
clusions from this study about transport of heat and mass in 
the bay are: 

a) The tidal flow is the dominant driving 
mechanism for general circulation in the 
bay. The tide flows primarily into and 
out of the southern part of the bay. 

The wind modifies the circulation to some 
extent though not significantly. 

b) The boundaries guide the tidal currents 
with velocities increasing as the bay 
narrows . 

c) The bottom topography affects the circula- 
tion by reducing velocities near the 
shallows. The Featherbed Bank significantly 
affects velocities around it. Water flow 
tends to go around the Featherbed Bank re- 
sulting in relatively high velocities in 
deeper regions adjacent to the banks . 

d) The tidal flow through the creeks produces 
only localized effects, 

e) The temperature distribution is predominantly 
determined by the bottom topography and is 
warmer in shallower regions. The Featherbed 
Banks shows up as a warm spot. Vertical 
diffusion is the dominant heat transfer 
mechanism. 

f) The effects of density change the temperature 
field over several hours. 
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g) The model has predicted within 1°C the 

3-D isotherms as compared to ground- truth 
data . 

It should be noted that these conclusions are 
for Biscayne Bay only, with its rather unique shallow and 
slowly varying depths, and not generalized conclusions for 
all basins. The aspect ratio, the sharpness of bottom slopes, 
and other factors, play a significant role. However, the 
model is general enough to accommodate a wide range of* geo- 
graphical and meteorological conditions. 
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A thermal plume discharging into a body of water 
is at first dominated by a jet like behavior. The jet then 
spreads with turbulent mixing at the edges. Depending on 
the densimetric Froude number, the bouyancy effect manifests 
itself by causing the plume to rise and spread at the surface. 
Beyond this region the jet spreads by the action of general 
"far field" circulation and diffusion. 

Fig. IV-50 shows the configuration of the surface 
discharge problem being formulated. The rigid-lid formulation 
described in previous sections still holds except that the 
boundary conditions are quite different. At the jet mouth 
the velocity, temperature and density profiles have to be 
specified. These are obtained from discharge measurements 
regularly recorded by power plant operating personnel. At 
the shoreline and bottom the same conditions as for far-field 
calculations apply, namely no slip, no normal velocity and 
adiabatic conditions. The outer boundary conditions are more 
difficult to specify. 

The domain outer boundaries are chosen such that 
the thermal anomalies caused by the plume are in the domain. 
The procedure for this is to first choose an outer boundary 
and then run the program; if no significant changes are 
observed at these boundaries the plume is considered enclosed 
in the domain. However, if the changes at the boundaries are 
significant, a larger domain is selected. The procedure is 
repeated until the thermal anomaly region is completely en- 
closed in the domain. 

At the outer boundaries the boundary conditions 
specified are zero gradients for velocity and temperature 
normal to the boundaries. This represents an approximation 
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that the boundaries are far enough so that gradients caused 
by the plume are non-significant. 

The method of solution is essentially the same 
as that for the far field. 

Preliminary results 

The Cutler Ridge plume has been approximated to a 
rectangular discharge 15 meters wide and 1.8 meters deep. 

The domain has a lateral extent of 300 meters and longitudinal 
extent of 600 meters as a first approximation. This is 
divided in 20 grid spacings in both horizontal coordinates. 

The vertical extent is divided into 4 spacings. The domain 
is assumed to be of a constant depth of 1.8 meters (same as 
discharge depth) . 

The discharge is taken to be at 25 cm/sec. The vertical 
and horizontal eddy viscosities are 1 cm /sec. and 1000 cm /sec. 
respectively. The basin is initially considered still. These 
conditions are chosen to test the program since the behavior 
of the jet is expected to be a symmetrical lateral spreading 
with the bottom boundary layer thickening V7ith axial distance. 

Fig, IV- 51 shows the surface' velocities for the dis- 
charge. The symmetrical spreading of the jet is clearly 
visible. The velocity decays away from the jet mouth as 
expected. 

Fig. IV- 52 shows the centerline velocity decay with axial 
distance. It can be seen that the decay is relatively small 
initially owing to inertial dominance of the jet. Then the 
lateral mixing causes faster decay. Finallyr the jet velocity 
asymtotically approaches zero since velocity gradients are 
small far av/ay from the jet mouth. 


The model results for the simple case considered 
indicate physically reasonable behavior. Effects of bottom 
topography, wind, temperature and currents are being in- 
corporated. The program already has the capability to in- 
clude effects of variable depth, wind, temperature and 
currents. Execution of more realistic cases has been 
initiated. 
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TABLE IV-I 


CASE 

WIND 
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HEAT TRANS- 
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BRUM 
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(°c) 

! 
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1 

10, SE 

1000 

26 

20 

3 

2 

10, SE 

500 

26 

20 

3 
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20, SE 

500 

26 
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10 , SE 

1000 

29 

20 

3 
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10, SE 
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29 

20 

1 
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10, SE 

750 

26 

20 
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10 , SE 

750 
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23 
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10, SE 

750 

30 
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Flow Chart for Free-Surface Model 
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Fig. IV-lO. Surface Velocities in Biscayne Bay With Wind From the North 
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Fig. IV- 11. Velocities at a Depth of 2 Meters With Wind From the North 
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Fig. IV- 12. Horizontal Grid System for Biscayne Bay 
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Pig. IV-14. Surface Velocities With Wind From the Southeast 
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Fig. IV-16. Velocities for Section J = 8 With Wind From the Southeast 
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Fig. IV- 17. Surface Currents with Incoming Tide 

in Biscayne Bay 
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Fig. lv-20. Currents at 2 Inters Depth During Tidal Flow Reversal in Bisoayne Bay 
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Fig. IV-21. Surface Currents for Ebb Tide in Biscayne Bay to Ocean 
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Fig. IV-22. Surface Velocities 
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Fig. IV-24. V^ocities at a Depth of 2 Meters 
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V. TURBIDITY STUDY 

A. Introduction 

Until very recently, the quality of our waters in 
lakes, streams and near the coast at sea has been deteriora- 
ting steadily. Dredging operations, discharge of industrial 
and municipal wastes, and agricultural usage all havve con- 
tributed to a decrease in quality of those very neighboring 
waters most useful to us for recreation and for fisheries. 

Since the general realization a few years ago that we must 
preserve our water heritage and limit our despoilation of this 
important natural resource, there seems to have been a dramatic 
turn of events. In some cases the rate of deterioration has 
apparently been slowed, and even reversed. Controls on 
activities which may serve adversely to affect water quality 
have multiplied, and considerable expenditures have been made 
or projected to improve water quality. Fishermen have re- 
ported improvements in the Great Lakes, and coastal waters 
such as Biscayne Bay seem somewhat clearer than a few years 
ago. But how is water quality to be measured? How are we to 
know whether we are winning the battle; whether the expenditures 
are fully effective. 

Of the measurements of water quality, that of turbidity 
or lack of clarity is probably the most useful. This is not 
only because measurements can be made optically, and therefore 
quickly without complicated chemical analysis, but because it 
is exactly the lack of transparency which turbidity implies 
that is perhaps the most direct and obvious measure of water 
deterioration. 

B . Defining Turbidity 

The literature discloses [1] at least nine distinct 
definitions of turbidity or lack of clarity. Here we shall 
consider that turbidity causes the cloudiness in water which 
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obscures objects or the bottom. We shall not be concerned 
with absorption, which limits the transparency of water but 
does not reduce the definition or contrast of what we see. 

The cloudiness that we wish to measure is caused by suspended 
small particles in the water. When light encounters these 
particles it is scattered in all directions. The result is 
a decrease in clarity and contrast of an image. 

Consider a beam of monochromatic light incident on 
a sample of water. It is assumed to have wavelength X and 
irradiance 1^. As it passes through the water, the intensity 
(irradiance) will decrease as energy is absorbed and scattered. 
As mentioned above, we will not concern ourselves with the ab- 
sorbed radiation, since absorption only decreases the intensity 
of the image, not its clarity or contrast. We will therefore 
evaluate only the scattering which occurs. 

In an infinitesimal volume d"v = dxdydz, a certain 
amount of light, dl, is scattered. Tnis amount is proportional 
to the incident intensity I^(X) , to the volume dxdydz, and is 
a function of angle and wavelength. 





Then 


dl = a(0/X) dx dy dz du) dX 


I =111/1 “(0-M 

X y z U5 X 


I dx dy dz duj dX 
o 


Here Qf(0,X) is the scattering coefficient which is to he 
measured. In those cases where the particles doing the 
scattering are small, of diameter d, the same order or smaller 
than the wavelength, X, of the light, the angular dependence 
is very complicated. In most of the turbidity which affects 
water quality, however, the particles are considerably larger 
than the wavelength, d » X, and the variations with angle 
are simple and not very strong; generally having the shape [2] 
shown in Figure V-2. 


20 40 60 80 100 120 140 160 180 

Scattering angle 9 

Figure V-2 
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Measurements of turbidity are thus equivalent to the 
measurement of scattering, which is dependent on angle as 
shown in Figure V-2. A type of scattering meter ^ the nephelo- 
meter is adapted to measurements at 90^. Other methods in- 
clude the taking of water samples, the Secci disk, and measure- 
ment of a light beam after passage through a definite length 
of water compared to its original intensity. These methods 
require slow sampling, or a bulky optical assembly, difficult 
to use in the sea, and liable to destruction by wave and ship 
motion. . 

C. Turbidity Study Plan 

We plan % three-stage approach to the problem of 
turbidity: 

1) The use of ERTS Satellite and/or airplane 
data suitably enhanced and applied as a time 
series for an overall view of a relatively 
large area . 

2) The use of a shipborne and later an airborne 
laser probe to measure backscatter turbidity 
as a function of depth. 

3) The development of a new ground truth turbidity 
meter based on a self-extinction method. 

These are discussed in turn: 

a) ERTS Satellite and/or Airplane Turbidity 
Measurements. 

Satellite and aerial photographic coverage aid in 
understanding the gross movements of suspended solids, providing 
information of circulation patterns [3], [4], [5]. When used in 
conjunction with field methods, such coverage minimized the 
physical sampling requirements. 

Perhaps more important than the photographic data 
taken by satellites and airplanes are the precise radiometric 
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data which are gathered by some satellites (ERTS-1 and 
ERTS-2) and specially equipped aircraft, Tiransmission of 

i 

this data are routed directly onto computer compatible tape. 

This increases the resolution by a factor of 7.9 over con- 
ventional photographic techniques. (Computer tape commonly 
has 7 bit resolution or 127 gray levels where photographic 
products have 16 gray levels. [6]). 

Bands 4 and 5 (approximately 5500 to 4500 and 6500 
to 5500 angstroms respectively) of the Earth Resources 
Technology Satellite-! (ERTS-1) have proved to be very effective 
in the observation of current patterns due to differential con- 
centrations of suspended sediment [7], Also observable by 
ERTS-1 imagery are waste effluents, production of planktonic 
organisms, river effluent and coastal landslides, variance 
in river discharge, and intensity of surf action [8], 

ERTS imagery photographic products may be enhanced 
using a densitometer, and the computer tapes enhanced using 
computer techniques. Enhancement techniques may be used to 
show the distribution of turbidity, from which a time series 
may be obtained from multiple imagery of the same location. 

Also the effects of tide may be observed but only to a limited 
degree since ERTS overpasses occur once every eighteen days. 

Using both enhancement techniques, photographic and 
computer, we propose to collect a time series of ERTS images, 
with corresponding field data, to draw conclusions into the 
nature, source and movement of turbidity through Biscayne Bay. 

By the use of time series, it is expectevd that the effects of 
turbidity may be separated from those due to reflection from 
the bottom, since the latter are probably much more slowly 
varying . 

b) Laser Application to Measure Vertical Sea 
Temperature and Turbidity 

This program is being developed by the Laboratory for 
Optics and Astrophysics of the Department of Physics, University 
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of Miami, under NASA contracts NASlO-8600 and NASlO-8795. 

It is described in the Detailed Technical Report 
NASA CR- 13 9 184 by J. G. Hirschberg, A. W. Wouters, F. N. 

Cooke, Jr., K. M. Simon and J, D. Byrne, dated January 1975. 

Tliis method evaluates the ratio of the backscattered 
Tyndall scattering to the backscattered Brillouin scattering. 
Only the Tyndall scattering causes turbidity, but since both 
are proportional to the incident light, the method is self- 
compensating . 

c) A Compact In-Situ Turbidity Measuring Meter 

In order to evaluate results from both preceding 
systems a compact ground truth turbidity measuring device is 
necessary. Ordinary available meters either possess a very 
small scattering volume and hence are insensitive, or are 
bulky and hazardous to operate at sea. 

This instrument measures the amount of backscatter, 
or scattering at about 180° per unit length of water. When 
looking downward into the sea with the light above, or when 
illuminating the sea with a light mounted close beside the 
viewer, this is the type of scattering which prevents observa- 
tion when present. In order to eliminate background inter- 
ference without a bulky mechanical target, near infrared 
radiation is used, produced by a solid state laser. This 
radiation is rapidly absorbed by water, so that the effect is 
like looking into a black hole at night, even in the daytime. 

The instrument contains one or more solid state lasers, 
emitting a collimated beam of radiation in the wavelength 
range from 700 to 850 nanometers. In this spectral region, 
the transparence of the water is quite low; from a meter to a 
few centimeters depending on the wavelength, see Figure V-3. 


) 


Absorption Coefficient in 



A wavelength is chosen for which the light has an 
effective path of 10 to 40 citia Thus, the background is black, 
exactly as though a black background target were used, without 
the necessity of providing a mechanical structure to hold it. 
This greatly simplified the apparatus for use at sea since no 
boom and disc are used, 

A further advantage is that the acceptance beam and 
illuminating beam are in the same volume, being anti-parallel 
so that the scattering volume is maximized. This makes the 
instrument far more sensitive than the nephelometer type where 
the beams are crossed. 
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Final ly, because the sea is absorbing to the 
radiation used beyond about 20 cm daylight obseir^rations are 
as favorable as nighttime for depths beyond about half a 
meter. 
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APPENDIX I 
EQUATION OF STATE: 

The generalized equation of state for water may be written as 
P = p (S,T,p) 
where p = density 

T = temperature 
P = pressure 
S = Salinity 

For most applications the effect of pressure on density is 
minimal since water is incompressible. Therefore , the significant 
variables are salinity and temperature. Thus i 
P = p (S,T) 

Dittmar (1884) showed that the majority of ions in sea water 
are present in remarkably constant ratios to one another. 

Salinity as defined by Forsch et al. (1902) is “the total 
solid material in grams contained in one kilogram of sea-water 
when all carbonate has been converted to oxide, the bromine and 
iodine replaced by chlorine, all organic matter completely 
oxidized, and the residue heated to a constant weight at 480°C." 

So, salinity is not the total salt content. Estimates by 
Lyman and Fleming (1940) related salinity S to total salt content 
by a linear relation : 

Stotai 0.043 + 1.0044 S. 

The chlorine-equivalent of the precipitated halides in grams 
per kilogram of sea water is called the chlorinity of sea-water. 

An empirical relation is : 

S =0.03 + 1.8050 Cl. 

here salinity S and chlorinity Cl are expressed in parts per 
thousand (0/00) . 
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Forsch et al, (1902) related specific gravity to 
chlorinity (at 0°C and atmospheric pressure) ; 

Oq = 0.069 + 1.4708 CL - 1.570 x 10”^Cl^ 

+ 3.98 X lO"^ Cl^ 

= specific gravity anomaly 

here = specific gravity. 

The effect of temperature variation is included as: 

2 

a. = A+Ba + CCT 
, too 

where 

A, B and C are functions of temperature. 

Knudsen (1901) used the formula : 

c-t = - So> + Bt<“o 

here 

2. = ~(0 -3.98)^(e + 283) 

^ 503.70(0 + 67.26) 

2q = -0.1324 

B = 1-A^. = 1-4.7867 x 1O"^0 + 9.8185 x 1O"^0^ 

-1.0843 X 1O”^0^ 

C = B^ = 1.8030 X 10“^6 - 8.164 x 10"^0^+ 1.667 x 1O"®0^ 

A = (4.531686 - 0.5459399^- 1.98248 x 10^^0^-1.438 x 10~^0^) 
(e +67.26) 

where 0 is temperature in degrees Celsius. 

For pure water - 2^ = 0.1324 
^ o o 

Therefore ^t " 

Another approach for particular cases is to fit the data 
to a suitable curve through density - temperature data. The 
C.R.C. Manual (1968) gives data for pure water in the range of 
10° - 50°C. - 
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A second degree least-square curve fit is used to obtain 
the empirical relation below: 

p = 1.000428 - .0000190 - .00000460^. 

Fig. IV-53 shows the variation of density with temperature 
and salinity calculated from the above formulae. 

In Biscayne Bay the salinity ranges from 35 to 38°/00. 

The corresponding density change has been ignored in the rigid-'? 
lid model since the rigid-lid model is primarily intended for 
inland water bodies where salinity is unimportant. Even in 
Biscayne Bay, the salinity variation is small so its effect on 
circulation is minimal and ignoring it in the thermal pollution 
math model can be justified. 

A second degree fit for the 35°/00 salinity data yields a 
formula relating density and temperature; which has been used 
in the rigid-lid model application to Biscayne Bay; 

p = 1.029431 -.0000200 -.00000480^. 

For the free surface model, the salinity may or may not 
be included depending on the densimetric Froude number of the 
water body. The system of equations described can then be 
used to approximate the effects of variations in salinity and 
temperature. 
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APPENDIX II 

A Model of the Sea Surface 
Consistent with Temperature Measurements 


Our measurements have shown that the Infrared measurements of the 
sea surface temperature are often not in agreement with measurements made 
with thermistors of the water just beneath the surface. In many cases 
these measurements are about a degree lower than the bulk measurements. 

When made from a satellite or high-flying aircraft, the reason for any 
discrepancy may possibly be ascribed to uncertainty in the quantity of 
water vapor in the atmospheric layer between the sensor and the sea sur- 
face. With our low-flying aircraft, NASA 6, this is not likely to be the 
case. Radiosonde data was available for all our experiments, and was 
utilized to correct the infrared observations. 

In addition, a hand-operated infrared sensor, the Wahl Corp. Heat < 
Spy was used directly from the boat. In this case where no water vapor 
correction applies, the measured sea surface was also usually about a degree 


cooler then the thermistor measurements.- In the following paragraphs a 
model yielding a possible explanation of these results is presented. 

As Figure 1 shows, the infrared in the 8 to 12 micron band which is 
used for infrared measurements of temperature is so strongly absorbed 
by water that infrared measurements are sensitive only to the upper 0.01 mm 
of the water depth. If for some reason this layer is warmer or cooler 
than the rest of the sea, the infrared results will not agree with bulk 
measurements. The below described model of the sea surface will serve to 
suggest why the top millimeter or so of the sea may be substantially cooler 
than the underlying water under typical conditions. We shall assume that 
these conditions include a relatively calm sunlit day, with the sun angle 
greater than about 30 degrees. The water surface is then relatively 

. " "" f 



smooth without foam, spray or many breaking waves. These conditions 
correspond rather closely to those obtaining during most of the experi- 
ments on which we report. 

The principal method of vertical heat transfer in the sea is by 
convection - the bulk transfer of fluid. However, unless the convective 
currents are very strong, or there is considerable wave motion, there 
will in general be a thin layer of water at the sea surface where con- 
vection cannot occur. This "skin" is maintained by the surface tension 
of the water, and heat is transfered through this layer by radiation and 
conduction alone. Let the thickness of the skin layer be t. In Figure 2, 
we represent the heat budget of the skin. A is the flux of sunlight 
falling on the sea. The sun's radiation is centered at about 530 micrometers 



and since water is relatively transparent to this radiation only a negligible 
part is absorbed in the skin; most of it passes into deeper water where it 
is finally absorbed or scattered back to the surface. B is infrared 
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re-radiated from the clouds and atmosphere, the total amount of which is 
small compared to A, except at night. C and D are the principal heat- 
loss processes for the sea. Let C be the infrared radiation, which for 
a temperature of 300°K has a maximum at about 12 micrometers. This cools 
directly only the upper 10 micrometers of the skin since water is 
practically opaque to this radiation as shown in Figure 1. D represents 
the heat loss by evaporation. This also cools only the upper layer of 
the skin directly, since only the few top molecular layers evaporate 
under the quiet sea-surface conditions we are considering. E is the 
delivery of heat from the deeper water to the skin layer, either by con- 
duction or convection. Finally F is the sunlight scattered by the water 
or reflected from the bottom. Typically this will be a small fraction of 
the arriving flux A. In this approximate treatment we shall neglect B 
and F since they are ordinarily, under the conditions assumed, much smaller 
than the other fluxes. 

For a steady state situation we can then write: 


and 


A = E 
E = C + D 


Now A is the arriving solar radiation. An approximate 24 hour average 
value for 25° latitude is^ 0.32 kilowatt/(meter)^. 

The principal mode of heat transfer in the skin layer is assumed to 
be conduction. The rate of heat conduction, dQ/dt, is given^ by 



- Ak grad T 


* By^nt. "Physical and Dynamical Meterology" Cambridge, England, 

1941, p. 155. 

^IngersoTl , Zobel and Ingersoll , "Heat Conduction" McGraw-Hill , 1948, p. 3 
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where A is the area, k the thermal conductivity, and T the temperature. 
In terms of the flux of heat, w, for conduction in one dimension this 
can be written: 



For water, k = 6 x 10'+ ergs/sec cm °IC‘+ , A ^ 0.32 kilowatt/m^ 
or 0.32 X 10^ ergs/cm^-sec 

£T _ 01 . 0.32 X 10^ 

3x " “ F 6 X lO** 

or, using differences: 

AT = 8ax . 

If AT is about 1®K, which is typical of our observations, ax is about 
1 millimeter, which seems reasonable, in view of the strong surface tension 
of seawater. This tends to explain the depression by about 1 degree K of 
the infrared temperatures , which are representative only of the upper 
1/100 mm of the ocean, compared to thermistor measurements at least a few 
millimeters below the surface. Refinements of this surface model may be 
made by combining it with the generalized force-surface model developed 
elsewhere in this report. 


‘^Quoted in Hirschberg, Wouters, Cooke, Simon and Byrne. Detailed 
Technical Report NASA CR-139184, January 1975, p. 32. 
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